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ABSTRACT 

In this paper a relationship between two approaches towards construction of genuinely 
two-dimensional upwind advection schemes is established. One of these approaches is of 
the control volume type applicable on structured cartesian meshes. It resulted (see [14], 
[15]) in the compact high resolution schemes capable of maintaining second order accuracy 
in both homogeneous and inhomogeneous cases. Another one is the fluctuation splitting 
approach (see [11], [3], [12], [17]), which is well suited for triangular (and possibly) un- 
structured meshes. Understanding the relationship between these two approaches allows us 
to formulate here a new fluctuation splitting high resolution (i.e. possible use of artificial 
compression, while maintaining positivity property) scheme. This scheme is shown to be 
linearity preserving in inhomogeneous as well as homogeneous cases. 
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1. Introduction 


The perfect scheme for numerical advection still awaits discovery. Its attributes should include 
high accuracy in regions of smooth flow, together with some positivity property that avoids over- 
and undershoots near discontinuities. It is well known that these two properties are not both 
attainable within the class of linear schemes. This implies that successful schemes must be non- 
linear, (i.e. data dependent). 

Thus the scheme must be furnished with a monitor function that measures the local smoothness 
of the data, for example by comparing consecutive gradients, as in the flux-limiter approach. This 
enlarges the stencil of the scheme beyond the minimum needed to attain a given formal accuracy. 
Such enlargement is undesirable for parallel implementation, and seems to impair the convergence 
of multigrid methods, as well as requiring supplementary conditions close to boundaries. 

In [15] a monitor function was introduced that compares gradients in different directions. This led 
to an advection scheme of the control-volume type with a minimally-enlarged stencil, and with good 
multigrid capabilities. Recently, there has also been developed the class of “fluctuation-splitting 
schemes. These schemes are intended to be coded as a loop over triangular (tetrahedral) elements, 
employing no data external to the triangle. Several non-linear variants of this technique have been 
devised that preserve positivity. 

In this paper we show a strong relationship between the two types of schemes. In fact, for linear 
advection on regular grids, certain schemes in each class turn out to be identical. This enables the 
transfer of insights and techniques from one class to the other. 

§2 presents some basic definitions and principles regarding the discretization of a conservation 
law. §3 documents several schemes of the control volume type, and shows how limiter functions 
can be introduced into them. §4 briefly describes the fluctuation-splitting approach, with some 
of the previously-employed positivity devices. §5 reformulates the nonlinear fluctuation-splitting 
schemes as limiter schemes and demonstrates their identity with some control-volume schemes. The 
performance of the fluctuation-splitting schemes is illustrated by some numerical experiments in 
§6. The conclusions of this work are presented in §7. 

The truly two-dimensional control-volume advection scheme that is capable of producing second 
order accurate steady-state solutions for the case of a non-zero source term is formulated follow- 
ing [15] in Appendix A. The high-resolution fluctuation-splitting counterpart of this scheme is 
constructed in Appendix B. 

2. Conservation law and its discretization 
C onsider the scalar conservation law 

(2.1) u t + V ■ f = s, 

where f — (/, g), and its nonconservative form, in which A = ( f u ,9u ) 

(2.2) u t + \ -Yu = s. 

In what follows, unless otherwise stated, we take A = (a, 5), with a, b constant 

(2.3) (d t + ad x + bd y )u = s. 

This assumption does not impair the generality of the construction. This is because using the 
conservative linearization procedure developed in [2] general nonlinear conservation laws can be 
represented locally by a linear constant coefficient advection equation of the form (2.3). 

We are interested here in the steady-state solutions (2.1), i.e. the case when d t u = 0. 

It will appear very useful to note that a steady-state solution of Eq.(2.1) in the homogeneous case 
(s = 0, the case considered through the entire paper, except for the Appendix) has the following 
property: it is constant along the direction of the convection velocity (characteristic direction). 
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Figure 1 . Control volume. 


The time stepping (forward Euler) procedure is considered here to be only a means to reach the 
steady-state. The solution update formula can be given by the following 

( 2 . 4 ) 

k 

i.e. the solution value in each grid point at the new time level can be represented as a combination 
of the solution values from the previous time level. 

Applying the TVD concept to characterize the schemes producing non-oscillatory solutions ap- 
pears to be too restrictive in two dimensions (see [4]). Therefore, we shall follow Spekreijse [16] 
and use the concept of positivity. 

Definition 2.1. A scheme is said to be of the positive type if the solution update formula can be 
written in the form (2.4) in such a way that c > 0, V&. 

Solutions obtained by the means of positive schemes will satisfy a certain maximum principle. 


3. Control volume approach 


Assume that the computational domain is divided into square cells (control volumes) each one 
associated with the cell-averaged value of the solution approximation (see Fig.l). Integrating 
Eq.(2.1) over a control volume and applying Gauss theorem we obtain 


(3.1) 



uX • dn = 



Since our purpose is to construct a second order accurate advection scheme, it is sufficient to 
approximate Eq.(3.1) using one- and two- dimensional mid-point quadrature rules. Thus, we obtain 


(^*2) t£ +1 — ^ [/+ — /- + 9+ *“ 9-] + SoA C 

where /+, f-,g + , are numerical fluxes computed at the centers of each side of the control volume 
Co. The question is how to compute them in order to obtain an advection scheme with desired 
properties. 

In order to simplify our further discussion in this section we assume that 


(3.3) 


0 < a < b 
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Figure 2. Stencils of the first order schemes: a) dimensional upwinding: b) N scheme. 


3.1. Dimensional upwinding. The simplest example of a positive linear scheme is given by the 
dimensionally upwinded scheme U with the following fluxes 


(3.4) 


fl = au3 

g^_ — bus 


This corresponds to the following update formula 

(3.5) < +1 = Uq - ~~ [6(uq - us) + a(uo - u 3 )] 

(Fig.2(a) presents the corresponding stencil). However, this scheme suffers from an excessive nu- 
merical diffusion. 


3.2. Zero cross-diffusion schemes. Here we present the zero cross-diffusion 2D scheme that was 
used in [14], [15] for constructing a nonlinear zero cross diffusion positive scheme (see also §3.4 and 
§3.5). The fluxes are 

,, f- D = au 3 - 2 6 (“3 - U 4 ) 

* ' ' gl D = bu 5 - l a(u 5 - u 4 ) 

These flux formulae can be motivated by considering the characteristics dy/dx = 6/a of (2.3). 
Specifically, produce the characteristic through the midpoint of the west face of the control volume, 
and find the value of u on it by linear interpolation (or extrapolation) in the interval 34. Multiply 
this by a to get /_. Similarly the characteristic through the midpoint of the south face carries a value 
found from the data 45, that gives g_ when multiplied by b. The additional terms found in these 
fluxes when compared with the U scheme can be regarded as antidiffusive terms which compensate 
for the diffusivity of the upwind scheme. It is natural to think of creating a high-resolution scheme 
by applying limiters to these terms. 

The update formula for this scheme is 

ii” +1 = u% - ^[^( u o - u 5 ) - ^(u 3 - u 4 ) + ^{u 5 - u 4 ) - ^(u 0 - u 3 )} 

(3.7) = u% - v, 4 - u 0 + us) + b{u°- u 5 ) + a (u B - u 4 )] 

(and its stencil is presented in Fig. 3(a). 

This scheme is not of the positive type - the coefficient multiplying u 3 in the update formula 
is negative in our representative case (3.3). Another defect of this scheme is that it will switch 
discontinuously when a or b changes sign. 
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A steady-state solution for the homogeneous advection equation obtained by means of a zero 
cross diffusion scheme will be second order accurate [14], [5]. 

It must be noted that there is an arbitrariness in the flux formulae. A term k(us — u 4 ) can 
be added to /_, and a term k(u 4 - u 5 ) to g_, without affecting the update formula, if k is any 
constant. In fact, one of the flux values is arbitrary. In the remainder of this paper, we adopt, as 
in [13], a convention that if b > a then g is always given by the characteristic interpolation (3.6), 
but if b < a it is / that is found in this way. 


3.3. N scheme* A positive, though first order accurate, modification of the scheme given by (3.6) 
is the following 

(3 8 ) f- - au 3 — 2 min (a, b)(u 3 - u 4 ) 

9 - = bus - \ min (a, b)(u 5 - u 4 ) 

or, because of our convention (3.3), 


(3 9) f- = au 3 - ±a(u 3 - u 4 ) 

9 - = bus ~ \a{u 5 - u 4 ) 

The corresponding update formula is (Fig.2(b) presents the corresponding stencil) 

(3-10) Uo +1 = Uq - ^-[b(u 0 - us) + a(u 5 - u 4 )\ 

This scheme was introduced by Rice and Schnipke [10] and was called the N (narrow) scheme in 
[14], [15] for its narrow stencil. Analysis in [13] shows that the N scheme is the optimal (in terms 
of cross diffusion) among the linear positive schemes in two dimensions. 

It is interesting to note that Raithby’s scheme [9] is in fact identical to the 2D scheme for the 
case b/2 < a < b. However, for 0 < a < 6/2 Raithby’s scheme would correspond to a blend of the 
2D scheme with the N scheme. The weight of the N scheme gradually increases with decreasing of 
a/6 and becomes equal to 1 when a = 0. 


3.4. S scheme. In order to combine the positivity property with second order accuracy a non- 
linear scheme has to be constructed. 

A positive zero cross diffusion scheme can be constructed in the way similar to TVD schemes 
in one dimensions. A “limited” zero cross diffusion correction is added to a first order accurate 
positive scheme. Such a scheme was developed in [14], [15] and was called the S scheme. It is a 
nonlinear modification of the scheme (3.6) and is given by the following fluxes 

(3.11) /- = f- ~ 2^(^543)(6 - a)(u 3 - u 4 ) 

9- = 9- 
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of a homogeneous advection equation (i.e. it has a zero cross diffusion) if is Lipschitz continuous 
and 


(3.18) ®(1) = h 

Remark 3.1 . The S scheme relies on a 5-point stencil (see Fig. 4(a)) - just one point more than the 
linear zero cross-diffusion scheme (3.6). This one extra point was needed to achieve positivity. 

It is remarkable that the choice of limiter function in order to obtain a positive scheme with 
zero cross-diffusion is dictated by relations (3.17) and (3.18) which are identical to those arising 
when constructing one-dimensional TVD-type schemes. Therefore, most of the limiters used in one 
dimension (see [18]) can be used in the present context as well. 

Remark 3.2. Another zero cross diffusion scheme was presented by Koren [7] 


(3.19) 

fl = /JV + I^( U5 _ U4 ) 

9- = 9- 

The update formula is 

u n+1 — u n - — 
“0 — “0 ^ 

2 -4^ (“5 «4) 

(3.20) 

+6(«o - U 5 ) 



+ 2 - «s) J 


The advantage of this scheme is that it switches continuously when a changes sign. 

It can be easily seen that this scheme is identical to the one-dimensional Lax-Wendroff scheme, 
assuming that y is the time-like direction. 

A detailed analysis of the class of zero cross diffusion schemes was done by Hirsch [5]. It is also 
proved there that a linear scheme with zero cross diffusion cannot be of the positive type, thus 
generalizing Godunov’s theorem to two dimensions. 


Remark 3.3 . An alternative route to a positive zero cross diffusion scheme is through the nonlinear 
modification of the scheme (3.19) 


(3.21) 



However, it is easy to see using (3.15) and assuming a symmetry property for the limiter, i.e. 


(3.22) 


*(A) 

R 




(and most of the commonly used limiters are symmetric, see [18]), that the scheme (3.21) is identical 
to the S scheme. 


Remark 3.4 ■ The choice of the ratio that the limiter function can rely on is non-unique. For 
instance, the upwind positive scheme with zero cross-diffusion presented by Hirsch k Van Ransbeeck 
[6] has a limiter function with the following ratio as its argument 


(3.23) 


#034 = 


ajup - M3) 
b(u 3 - u 4 ) 
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3.5. S* scheme. One of the main objectives of this work is to establish links between the control 
volume and fluctuation-splitting approaches towards the construction of the truly two-dimensional 
advection schemes. It is important for this purpose to introduce the scheme S* which has a more 
‘triangular flavor’. The limiter schemes above are based on ratios R{ 3 k of gradients in orthogonal 
directions. Now we introduce the ratio 7?* jfc , defined as 


(3.24) 


-a(ttt - uj) 
ijk (b-a)(u k -u : ) 


where the pair of points j,k are adjacent along a diagonal. The fraction has been normalized so 
that we still have R* jk ~ 1 for smooth, nearly steady solutions. 

The S* scheme can be defined by the following numerical fluxes 


(3.25) 

with 

(3.26) 


f- - f--\V(Ro 4 3 ){b - a){u 3 - u 4 ) 


9 - = 9- 


R 


a(uo — uf) 


043 


:)• 


(6 - a)(u 3 - u 4 ) 

The update formula in this case is (the stencil is presented in Fig.4(b)) 

b(u 0 - -us) + «(«5 - “ 4 ) 

At I 

(3.27) 


n + 1 


U n - 


+\^(R 043 ){b - a)(u 3 - u 4 ) 

L -^(/?75o)(h - a)(«o - “ 5 )] J 


Note that the following identity holds 

(3.28) 

Lemma 3.5. If the limiter 'I' = R ) satisfies the following inequality 

(3.29) - _ R 
then the S* scheme is of positive type . 


*(Ro 43 )(b - a)(u 3 - u 4 ) = - ^§^- a(u 0 - u 4 ). 

* 7(343 


0 < < 2,0 < 4 >(R) < 2. 


Proof. Using the identity (3.28) the update formula for the S* scheme can be rewritten as follows 


(3.30) 


u ; 


n+l 


Un - 


At 

~h 


( b — a) (l - ^ 4* (T? 75 o)) (“o - U 5 ) 

+« 0 - (“° " J 


□ 


It is obvious from (3.30) that the scheme is positive if the inequality (3.29) holds. 

Lemma 3.6. 7/^(7?) is Lipschitz continuous and 
(3.31) *(1) = T 

then the S* scheme has a zero cross diffusion. 

Proof. This lemma is a simple corollary of Lemma A.l concerning the more general scheme con- 
structed for the inhomogeneous equation. □ 
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Figure 5. Illustration for construction of fluctuation-splitting schemes. 


4. Fluctuation splitting approach 

In this section we present a brief description of a fluctuation splitting approach. For more details 
see [3], [12], [17]. The schemes are designed for use on unstructured grids and their description is 
initially given in that context. Later, to compare with the control- volume approach, we specialize 
to structured grids, and later still we reconsider the general case. 

Consider a numerical solution of the two dimensional linear constant coefficient advection equa- 
tion 


(4.1) 


u t + A • Vu — 5. 


The grid is taken to be an arbitrary triangulation of the domain. A typical element T of such a 
grid is shown on Fig. 5. The integral of u t over the triangle T is called the “fluctuation” 


(4.2) 


'-JL 


Ufdxdy . 


Taking into account Eq.(4.1) and applying Gauss’ theorem we obtain 


(4.3) 




(p uX • dn — I I s dxdy , 
JdT J JT 


?dT J JT 

where dT is the boundary of the triangle T and dn is the inward scaled normal to an element of 
the boundary. Assuming that the source term s varies linearly within triangle T we get 

(4.4) a T = JJt s dxdy = Y] s i , 

where Sj is the area of the triangle T. We shall return to the inhomogeneous case in Appendix B. 
Here we assume that s ~ 0 (and therefore aj = 0). 

(4.5) u t + A • Vu = 0. 

Assuming also that u varies linearly within triangle T we get 


(4.6) 


T L — * 1 — » 1 — » 

4 > — T ^2)^ * + “( u 2 + ^3)^ * + “(tti + W3) A * 7?2, 
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Figure 6. Median dual cell 


where 

(4.7) 



and n t is the inward normal to the side 2?, scaled with its length. Note that 

3 

( 4 . 8 ) = 0 

i = l 


and therefore 
(4.9) 


3 


J> = 0. 


This allows (4.6) to be rearranged as 
(4.10) 4> t 
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i-l 


Several alternative expressions can be obtained for using (4.9), for example, 

<t> T = -k 2 {u 2 - ui) - k 3 (u 3 - ui) 

— -k 3 (u 3 - u 2 ) - ki(ui - u 2 ) 

(4.11) = -fci(tti - u 3 ) - k 2 (u 2 - u 3 ) 

The computed fluctuation should then be distributed (split) among the vertices of the triangle T 

(4.12) := SiUi + ajAt4> T , i = 1,2,3, 

where S* is one-third the area of all the triangles having V, as a vertex (the area of so-called median 
dual cell of area 5, around i, see Fig. 6) and 


(4.13) 


3 



The latter is necessary for the discretization to be conservative (see [3]). 

After adding contributions from all the triangles T having a common vertex at point i , we obtain 
the following scheme 


(4.14) 


i, 71 + 1 


= K + ^7 X] °?4 >t 


Si 
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Figure 7. Two possible triangle orientation with respect to the flow direction. 


We assume here that each triangle sends contributions to its own vertices only. Therefore, the 
resulting scheme for u t has a compact stencil restricted to at most the immediate neighbours of m. 
The remaining question is how to chose af. 

An additional consideration is positivity of the constructed scheme. We will define a more 
restrictive condition of local positivity that requires contributions of each triangle separately to 
be positive. 

Definition 4.1. In equation(4.12) choose from the preceding list the i th definition of (j> T , so that 

(4-15) S t u t := SiUi - ai[kj(uj - u t ) -f kk(uk - «,•)]. 

Then the scheme is said to be locally positive if the three quantities 

a t kjAt, otik^At, 1 -p o^^kj T kk)At 

are all positive 

This condition is easier to implement because it retains the basic property of a fluctuation- 
splitting scheme that all triangles can be processed without reference to any other data. However, 
it is more restrictive than necessary because it does not recognise that compensating changes from 
other triangles may restore positivity. It is the condition currently used in design of fluctuation 
splitting schemes. Application of the non-local positivity condition may allow addition of artificial 
compression to the scheme. This will be discussed in §5.2. 

Definition 4.2. The fluctuation-splitting scheme is called linearity preserving if whenever the fluc- 
tuation on the triangle T vanishes then the scheme leads to a zero update in each of the three vertices 
of the triangle. 

It was observed by Deconinck et al [3] that linearity-preserving schemes give solutions that are 
second-order accurate in the steady state. 

4.1. N scheme. This linear, positive, fluctuation splitting scheme is called the N scheme because 
it was found by Roe [12] that it is closely related to the control volume N scheme. (We shall return 
to this point in §5.1.1). Note that there exist two possibilities regarding ki for i = 1,2,3 (because 
of Eq.(4.9)): either one of them is positive and the rest are negative or two of them are positive 
and one negative. These two situations correspond to a triangle with either one inflow face or two 
inflow faces. Without loss of generality we can consider the following two cases (see Fig. 7) 

(a) : ki < 0, k 2 < 0, k 3 > 0, 

(b) : ki > 0, k 2 > 0, k 3 < 0. 
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Figure 8. Construction of the NN scheme. 


The obvious choice for the first case (Type I triangle) is to send the entire fluctuation to the 
downstream vertex 3 (i.e. setting aj = — 0, aj = 1 ) 

(4.16) S 3 U 3 +* = 53^3 “ — U 3 ) + £ 2(^2 — ^ 3 )] 

Although interesting schemes can be created that do not follow this rule, we will assume for the 
remainder of this paper that all Type I triangles are dealt with like this. 

Case (b) is more complicated. If the fluctuation is distributed according to the following formulae 

Siu? +1 = Siu? - At[&i(ti! - u 3 )] 

^ ' ' S 2 u n 2 +' = S 2 v% - At[k 2 (u 2 - u 3 )] 

the resulting scheme is obviously positive and this will define the N sceme. Strong convergence of 
the N scheme on general triangulations is proved by Perthame [8]. However, it does not have the 
LP property. This is because the contributions to the residuals at both vertices 1 and 2 may be 
different from zero even if the fluctuation on triangle T vanishes. 

4.2. NN scheme. It was pointed out in [11], [3], [12] that replacing the advection velocity A by 
its normal to the level lines component A n normal to the level lines (or tangential to the gradient 
direction) of u within triangle T 

- w Vu 

(4.18) A n = (A-Vu)--^— 

I \u\ l 

does not affect the residual. It is also obvious that substituting the advection velocity by the 
following 

(4.19) A* = A n T OX 
does not affect the residual either. Defining 

, * lr* . 

(4.20) k* = - A • n t 
and 

5x< +1 =s 1 ti?-At[*;(tti-« 3 )] 

^ ’ 52«2 +1 — S 2 u^ - At[kl(u 2 - u^,)) 

one can observe that if 0 — 0, the NN scheme (4.21) will be linearity preserving. This is because 
A n vanishes at the steady state (and therefore fcj and k% vanish in this case). However, vector A n 
may extend out of triangle T (Fig. 8). Therefore, the scheme in this case is not positive, i.e. one of 

ii 



the s in (4.21) will be negative. To obtain positivity of the scheme in this case 6 should be taken 
such that the vector A* will be directed along the face of triangle T. This means that the entire 
residual should be sent to only one of the outflow nodes. 

Because the effective advection speed is data dependent (when 6 = 0 it is the drift velocity of 
the contours) the NN scheme is nonlinear in a way that permits it to escape Godunov’s Theorem. 
Although the mechanism involved appears very different from that employed by limiter schemes, 
we will see in the next section that a unification is possible. 

It was pointed out in [3] that when the solution is almost a constant, the expression for X n (4.18) 
becomes ill-defined, which may affect the convergence. A way to overcome this difficulty suggested 
in [3] was to switch gradually from A* to A in regions where the gradients are small. 


4.3. Level scheme. It was argued by Roe in [12] that the NN scheme is not optimal in terms of 
the maximal allowed time step and another scheme for Type II triangles was suggested there. This 
scheme was later named the Level scheme. 

A brief description of the algorithm is given here. The reader is referred to [12] for the complete 
derivation. 

The residual is distributed among the nodes according to the following formulae: 

U 22) Siu^ +1 = + A ta^ 

1 = S 2 U 2 + A.tQ'2<f>T 

It is very interesting for the discussion in §5.1.3 to note that in [12] the following quantity 


(4.23) 


A 


^1 ~ «3 

u 2 - U 3 


was suggested to distinguish between two cases 

(a) : A > 0, i.e. the level line passes through triangle T, 

(b) : A < 0, i.e. the level line lies outside of triangle T. 

In case (a) only one vertex is updated: the vertex that the characteristic line passes closer to it 
than the level line. 

In case (b) both vertices are updated. The residual is distributed according to the following 
formulae 


ai = Si(m-tt 3 ) 

(4 24 ) Si(u 1 -'U3)+S 2 (u2-u 3 ) 

' ‘ ' ry = S 2 (u 2 -u 3 ) 

2 Si (tzi - u 3 ) + S 2 (112-^3) 

This choice of ct’s is argued in [12] to be optimal in terms of the maximal allowed time step. Also, 
this choice preserves the direction of the level lines; hence the name of the scheme. 

However, Deconinck [1] has observed that the Level scheme may switch discontinuously if A is 
almost parallel to one side of the triangle, say side 1, but (u 2 ~u 3 )/(ui -u 3 ) > 0. This demonstrates 
that A was an unhappy choice of parameter. 


5. Unified approach 

It is shown in this section that there is more in common between the two previously reviewed 
approaches than just a philosophy to make constructive use of nonlinearity. In fact the fluctuation 
splitting methodology and control volume compact scheme methodologies can be seen as two dual 
ways to interpret the same approach. 

Consider the following scheme on Type II triangles 

‘Sl? i i l+1 = S\U y + A£[/ji(tti — U3) + ^ {R\32)k 2 {u 2 — Us)] 

S2 U-2 + ~ S2U2 + At[k 2 (u 2 — Us) — ^ (R\ 32 )k 2 {u 2 — W3)] 
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where Rij k now compares contributions to the fluctuation from two sides of the triangle, i.e. 

ki(ui — Uj) 

k k (u k -uj)' 

Again, we have R tjk — 1 for smooth, nearly steady flow. 

kjjui ~ u 3 ) 
k 2 {u 2 - u 3 ) 

The conservation property of this scheme is obvious from (5.1). 

We have the identity 

_ ^(# 132 ) 


(5.2) 


(5.3) 


Rijk — 


#132 — ~ ' 


(5.4) 


*(#132)M«2 “ U 3 ) = — 


R 


132 


'M U 1 - U 3 ) 


Theorem 5.1. If the limiter 'fl(P) satisfies the following inequality 

m) 


(5.5) 


0 < 


R 


< 1, 0 < V(R) < 1, 


(5.6) 


then the NNL scheme is (locally) positive. 

Proof. Using the identity (5.4) the scheme (5.1) can be rewritten in the following form 

Siu" +1 = + At[l - *^3^ 3M«i ~ u 3 ) 

S 2 u” +1 = S 2 m£ + At[l - ^ (# 132 ) 1 ^ 2(^2 - “ 3 ) 

It can be concluded from the latter that the NNL scheme is positive if inequality (5.5) holds. □ 

Remark 5.2. By contrast to the control volume schemes, it is here possible only to use limiters for 
which the upper bound on 'I','!'/#, is equal to 1.0 rather than 2.0. That is, the limiter cannot 
be compressive. We show later that compressive limiters can be allowed within the fluctuation 
splitting scheme given sufficient information about mesh connectivity. 

Theorem 5.3. The NNL scheme is linearity preserving if 

(5.7) |#(#)-1|<C|#-1|. 

for some positive constant C . 

Proof. We have 

Siu ” +1 — Sin" = 5u\ = Af[fci(tii — 113 ) + *(Ri 32 )k 2 (u 2 — 113 )] = P + 'L(#)Q 


(5.8) 


ibgi 


S 2 u 2 +1 - S 2 M 2 = $ u 2 = At[k 2 {u 2 - u 3 ) + Ki32 
We wish to show that if P + Q = 0 then 6u j = 8u 2 = 0, or equivalently 

|<5u! - Sufi < K\5u\ + Sufi, 
for some positive constant K. Now, if 

I* - 1| < C\R- 1| 

2|tf-l||Q| < 2C|P + Q|, 

|P + Q| + 2|¥-l||Qi < (2C+1)|P + Q|, 
|P + Q + 2(#- 1)Q| < (2C + 1)|P + Q\, 
|P-Q + 2*Q| < (2C+1)|P + Q|. 


^ 1 (m 1 -m 3 )] =Q-y(R)Q 


then 


(5.9) 


This last statement translates, with K = 2C + 1, into the one we want to prove. 
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2 1 8 



Figure 9. Possible triangulation of the Cartesian grid. 


2 1 8 



4 5 6 

Figure 10. Possible triangulation of the Cartesian grid. 

Remark 54 . It is interesting to note that the NNL scheme with the “minmod” limiter is identical 
to the Level scheme in the case when the level line lies inside the triangle. To illustrate this assume 
that Rim > 1. For the minmod limiter this means that ^(/?i 32 ) = 1. It follows from (5.1) that the 
entire fluctuation contributes to the solution update at the node 1 - exactly what the Level scheme 
corresponds to. 

5.1. Links to the control volume approach. 

5.1.1. N schemes. Consider the grid as illustrated in Fig.9. The fluctuations for the triangles 
A,B,C are given by the following expressions 

4> A = -f[(£> - a)(u 3 - u 4 ) + a(iio- u 4 ) ) 

4> b = — f [a (up - n 4 ) + (6- a) ( M q - n 5 )] 

4> C = [ a («7 - tt5) + (fr-q)(tto-t* 5 ) ] 

The underlined terms here contribute to the residual at the node 0 according to (4. 16), (4. 17). 
Assembling all the contributions the following update is obtained 

(5*10) < +1 = + —[b(u$ — uq) + a(u 4 — uq)] 

It was pointed out by Roe [12] that (5.10) is identical in this case to the update formula for the 
control volume type N scheme (3.10). That is why these two schemes are called N schemes. Another 
interesting observation made in [12] is that if the N scheme is constructed on the triangulation 
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illustrated on Fig. 10 


4> A = -|[(6 - a)(u 3 - u A ) + ajup - u 4 ) ] 

(f> B = - 1 [ a(tt 0 - u 4 ) + ( b - q)(uq - tr 5 ) ] 

4>' 0 - -|[o( u 7 - u 5 ) + (6 - a)(Mp - U 5 ) ] 

it can be seen from the update formula 

(5.11) u£ +1 = u o + ~ u ° ) + a ( U3 ~ v o)] 

that the scheme in this case is identical to the dimensional upwinding (3.5). 


5.1.2. 2D scheme. Consider again triangulation illustrated on Fig. 9. The fluctuations on the tri- 
angles A, B,C are again given as follows 

<P A = -| [(j> - a)(«3 - ttj) + a(qo - q 4 ) ] 

<t> B = — ? [a(ttp - Uj) + (6 - q)(uo - -tt 5 ) ] 

<tf = -%[a(u7 - u 5 ) + (b - a)(u 0 - U 5 )] 

Note, that now only triangles A. B participate in the molecule at node 0, however they contribute 
their entire fluctuations. This would correspond, in the fluctuation-splitting approach, to a rule 
that each triangle only updates it’s most downwind node (the one for which A • x is maximum). 
The update formula in this case 

, . _ A t r (a — b), . (b - a) , (a + b) 

(5.12) U 0 + = U 0 + — [ (u 3 - Uo) -I («5 - u o) H ^ (lt 4 — T* 0 )] 

is identical to the 2D control volume scheme (3.7). 


5.1.3. Linearity preserving positive schemes. Consider triangles A,B,C as illustrated in Fig. 9. De- 
fine 


(5.13) 


7?043 — 


a(uo - u 4 ) 

( b - a)(u 3 - u A )' 


R 750 — 


a(uj - u 5 ) 

(6 - a)(u 0 - u 5 ) 


The fluctuations on these triangles are 


4> a = -|{[1 - ^(/?043)](^ - a)(u 3 - U 4 ) + ajup - m 4 ) + ^(# 043 ) (6 ~ a)(g3 ~ ^ 4 ) } 
(5.14) <f> B = — |{ a(u 0 — u 4 ) + (6 - a)(ttp - us) } 

4> c = -|{o(ti 7 - U 5 ) + ^(/?75o)(^ - a){uo - U 5 ) -I- [ 1 - ^(fi75p)](ft ~ < 0(^0 ~ ** 5 ) } 


The underlined terms contribute to the residual at the point 0 

a, [ {{b-a)-^(R 75 o)(b-a))(u 5 -u 0 ) 
(5.15) «S +1=u o + X + (a+ i'V{Ro43){b- a)) (u 4 - u 0 ) 

[ +{-l'H(Ro43){b-a)) (u 3 -u 0 ) 

which is identical to the update formula for the S* scheme (3.27)! 


5.2. Artificial compression. Note that it was possible to use a compressive limiter 2 * for the 
S* scheme. This was possible to show by requiring non-local positivity, which is a less restrictive 
requirement. If we were to restrict our attention only to structured grids, we could use the equiva- 
lence of the S* and N N L schemes to justify using compressive limiters in the NNL scheme also. 
We now 7 discuss how compressive limiters may be employed w 7 hen we know 7 something about the 
connectivity of the grid. 

Consider triangle C (Fig. 9). Collect from triangles B and C those terms in the fluctuations that 
contribute to the update at node 0 and contain a factor (us — uq). 


2 A compressive limiter can be defined as one for which the upper limit on $ / R is greater than unity. 
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Figure 11, General triangulation. 

^o+<£o = ~ ^(^75°)](& ~ a)(«5 - «o) + (b - a)(u 5 - u 0 )} 

(5.16) = ~ [2 — 4'(i?75o)](6 — a)(w 5 — uo) 

Similarly, collect from triangles D and C those terms in the fluctuations that contribute to the 
update at node 7 and contain a factor (^5 — uj ). 

<t>7 +<t>7 = ~ a )( u 5 - u 7 ) + (6 - a)(u 5 - u 7 )} 

(5.17) = ^ [2 _£&2) ](6 _ a)(M5 _ U7) 

^ J1750 

It is obvious that both these contributions will be positive if the following inequality holds 

(5.18) 0 < <2, 0 < V{R) < 2. 

This means that a compressive limiter can be used in this case. However, this case is particularly 
simple: the grid is nicely structured and the advection velocity is constant. The case of unstructured 
grids and/or variable advection velocity is more complicated. ( The nonlinear case can be reduced 
for this purpose to the variable velocity case by using the multidimensional conservative linearization 
(see [2])). 

Consider this general situation as illustrated in Fig.ll. Assume that triangle C is of Type II, 
while the neighboring (through the inflow faces of C) triangles B and C are of Type I. Collect, as 
before, the contributions to the updates at nodes 0 and 7 containing factors (u 5 — uq) and (n 5 — u 7 ). 

4>o + 4>o = ^{[1 _ ^(^so)]^^ — u o) + &<f(M 5 — ^o)} 

(5-19) = + k§) - 5o)A£]K - u 0 ) 

4>7 +<f>7 = - U 7 ) + kj («5 - U 7 )} 

( 5 - 20 ) = §[(*? + 

z ^750 
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Figure 12. Positivity region for the limiter used in the NNL scheme. 


where k ’ s are defined according to (4.7) (superscripts indicate the triangle). Defining the following 
Artificial compression factor' 


(5.21) 


. k% + kf + kfi 

7ae = min(- — , c — 

k 7 *0 


) 


we can state 


Theorem 5.5. The NNL scheme is of the positive type if the limiter ty(R) employed in triangle C 
satisfies the following inequality 

(5.22) 0 < < Tac, 0 < V(R) < 7a c- 

Proof. The proof is obvious from the previous argument. □ 

The positivity region (5.22) for the limiter function used with the NNL scheme is illustrated on 
Fig. 12. A class of limiters defined by 

(5.23) V* = max(0, min (4>R, 1), min (R, <f>)) 
with the parameter 4> satisfying 

(5.24) 1 < <P < lac- 
appears to suit this case very well. It is obvious that the choice 

(5.25) f> = 1 

gives the minmod limiter which corresponds to the lower bound of the positivity region (Fig. 12). 
The choice 


(5.26) 4> = 7ac, 

results in the “generalized” highly compressive Roe superbee limiter and corresponds to the upper 
bound of the positivity region (see Fig. 12). 

The amount of artificial compression that the scheme is allowed to use should probably be 
problem dependent. Usually it does not make much difference for the resolution of strong shocks 
and may lead to some undesirable effects in the smooth regions (for example, an advected sinusoidal 
wave may start looking like a square wave). However, it can make a great difference in resolving 
contact discontinuities. 
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Figure 13. Regular upwinding and N scheme. 

6. Some numerical experiments 

In this section we present some numerical experiments which illustrate the performance of the 
genuinely two-dimensional advection schemes. Two important aspects are considered: resolution of 
discontinuities and the accuracy in smooth regions. All testcases presented in this section consider 
the linear advection equation 


Ut -\- Uy - 1- .5 u x — 0. 

6.1. Resolution of discontinuities* An extensive set of numerical experiments with the S scheme 
concerning the resolution of discontinuities is presented in [15]. Here we shall illustrate this ability 
of the S* (or NNL) in comparison to the 5 scheme. 

We consider here the advection equation (6) on the square [0, 1] X [0, 1]. The boundary conditions 
and the exact solution steady-state solution are given by 

u = H(y- (2x - 2)), 

where H is the Heaviside function. The grid used is 50 x 50 points. 

Figuresl3-16 present contour plots of the numerical solution to this problem obtained using 
different schemes. The performance of the first order regular upwind and the N schemes is illustrated 
by Fig. 13. It can be observed that the N scheme is superior to the dimensional upwinding in terms 
of discontinuity resolution. However, the width of the numerical layer representing the discontinuity 
is 0{h 1 ! 2 ) in both cases. 

Fig. 14 presents the results obtained by the S and the S* (which is the same as NNL scheme in 
this case) employing the non-compressive minmod limiter. Fig.15 corresponds to the same schemes 
except that the van Leer (harmonic mean) limiter is used. Experiments with the highly compressive 
superbee limiter are presented in Fig. 16. It can be seen that the discontinuities are represented by 
sharp layers with no oscillations in the numerical obtained by both the S and the S* schemes with 
all the three different limiters. Another conclusion is that the artificial compression may have a 
strong effect on the sharpness of such a layer. It can also be observed that the S* is slightly superior 
to the S scheme in resolving contact discontinuities. This can be explained by its slightly narrower 
stencil. 


18 





Figure 14. S and NNL scheme with “minmod” limiter. 



Figure 15. S and NNL scheme with van Leer limiter. 

6.2. Accuracy in smooth regions. The purpose of the experiments presented here is to verify the 
accuracy of the S* (or the NNL scheme on structured grid) in the smooth regions and to compare 
it to that of the S scheme. An extensive set of numerical experiments with the 5 employing 
different limiters was presented in [15]. Here we restrict our attention to the minmod limiter since 
the S* scheme (as well as the NNL scheme in general unstructured grid case) satisfies the local 
positivity property in this case. Therefore, the treatment of each triangle is purely local. This 
results in a significantly more economical numerical algorithm, though at the expense of slightly 
inferior resolution of contact discontinuities. 

The testcase considered is again the advection equation (6) on the square [0, 1] X [0, 1]. The 
boundary conditions (and the exact steady-state solution) are given by 

u = sin(4:r - 2 y) 
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Figure 16. S and NNL scheme with Roe superbee limiter. 


S scheme 
1.40 x 10~ : 
3.38 x lO"' 
8.54 x IQ- 4 


Table 1 . L i solution error norm (measured on the [0,0.75] x [0,0.75] subdomain to 
avoid the influence of numerical boundary layers) for the problem u t + u y + .5 u x = 0 
with the boundary conditions (and the exact steady-state solution) u = sin (4a; -2y). 


The Li solution error norm measured on the subdomain [0,0.75] x [0,0.75] (to avoid any influence 
of the numerical boundary layers) is presented on Table 1 for the cases of the 5 and 5* schemes 
and three different levels of resolution. 

The experiments with both schemes demonstrate second order accuracy at the steady-state of 
the both schemes. The solution error in the case of the S* scheme is slightly smaller than for the 
S scheme. The explanation for this is the slightly narrower stencil of the 5* scheme. 

7. Discussion and conclusions 

This work establishes a link between two types of genuinely two-dimensional advection schemes. 
The truly two-dimensional schemes of the control volume type were constructed in [14], [15], The 
fluctuation-splitting type schemes were introduced in [11], [3], [12], [17]. The 5* scheme presented 
in this work can be interpreted as a representative of both classes. This allowed us to transfer 
ideas and insights between the two approaches and to combine the desirable properties of both 
classes in a single scheme. This is accomplished by the NNL scheme, which on one hand (as a 
fluctuation-splitting scheme) is formulated on unstructured triangular grids. Thus, it allows for a 
simple treatment of complex geometries. On the other hand the NNL scheme inherits the following 
properties of the control volume type genuinely multidimensional advection schemes: 

• to use regular limiter functions in order to formulate the fluctuation-splitting schemes for 
unstructured triangular meshes; 

• possibility to use artificial compression, which can improve the resolution of discontinuities: 
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• to formulate a fluctuation-splitting scheme having the LP property for the advection equation 
with nonzero forcing term. 

One of the important objectives of this research direction is to construct highly efficient and 
robust multigrid steady-state solver for the compressible flow problems. The main motivation for 
constructing the multidimensional schemes of the control volume type in [14], [15] was to obtain a 
high-resolution scheme such that Gauss-Seidel relaxation (the simplest and very efficient smoother) 
is stable when applied directly to to the resulting discrete equations (which is not the case for 
dimensionally-split high-resolution schemes [16]). This made it possible to construct a highly effi- 
cient multigrid steady-state solver. The genuinely two-dimensional advection schemes formulated 
in this work ( S * and N N L) inherit this property from their control volume type predecessors. 
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Appendix. 


Appendix A. S scheme - general inhomogeneous case 


Consider the following scheme 


(A.l) 



aus — o(k (^3 - ^ 4 ) — hs{_) 

bus — - U4) — hs 9 _ ). 


It has been shown in [14], [15] that the 2D scheme is second order accurate at steady state for 
inhomogeneous equation. The definition of s/_,s 9 _ is non-unique (see [14], [15]). We consider here 
one particular choice which appears to be of particular importance for the purpose of this work 
(see Appendix B): 


(A.2) 


f 2 1 

S -~3 S3+ 3 S4 


, 2 1 

(A. 3 ) s_ — -S5 + -S4 

If one is content to obtain first order accuracy the previously defined N scheme can be used. 
However, the purpose of it here is to serve as an intermediate step towards the construction of the 
second order accurate S* scheme. Therefore, we shall modify it. For our representative case a < b, 
putting 

(A.4) /3 = min(l,^) = | < 1, 


s~ — ~(so + s 3 + S-l) 

au 3 - j{a(u 3 - u 4 ) - h[s f _ - (1 - /?)«-]} 
bu 5 — |(a(tt 5 — u 4 ) — hs 3 _). 


and defining 
(A.5) 

the N scheme can be given by: 

f N _ 

(A-6) A' _ 


The S* scheme is defined by 


(A.7) 



/ s _ - ly{Ro43)((l> ~ a ){u 3 - u 4 ) - (1 - ( 3 )hs ) 
</ S - 


where 

(A. 8 ) R 043 

The update formula in this case 

« 5 +1 = 


(A.9) 


a(u 0 — U4) — ( 3 hs _ 

( b — a)(u 3 — ^4) - (1 — j 3 )hs _ 

b(u$ - uq) + a(u 4 - u 3 ) 

-l'J'(fl 7 5o)[(6 - a)(u 5 - Uq) + h{ 1 - /?)s_] 

+-^(R 0 43)[(b - a)(u 4 - u 3 ) + h( 1 - /?)«+] 

+^[s+ - st + - s 9 + + (1 - / 3 )(s+ - s_)]} 
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or 


u o +X ~ u o + ~j^{ b{u 5 - u 0 ) + a(u 4 - u 5 ) 

— (i?75o)[(6 — a)(u 5 — u 0 ) + h( 1 - /3)s_] 
+-^(i ?° 43 )[(6 — a)(u 4 - U 3 ) + A(1 — /?)s + ] 

(A-10) + ^h [2s ° + (1 - + (2 - 0)*4 + (1 + P)s 5 + /fc 7 ]} 

Lemma A.l. If ^ — \l/(i?) is Lipschitz continuous and 
(A- 11) *(1) = 1, 

then the S* scheme is second order accurate. 

Proof. Note that 

(/+ - /-) - ul D - f- D ) 

= (1 - V(R75o)){(b- a)(u 0 - u 5 ) - (3hs+) 

— (1 - ^(/?043))((& — ffl)(«3 — V-4) - /3hs _) 

= h(V(Ro 43 ) - V{R 750 ))({b - a)(u y - .5 hu yy ) - (3s) 

+h 2 { 1 - ^(^043))((6 - a)u xy - (3s x ) 

(A. 12) +0{h 3 ) 

The second order accuracy of the S* scheme requires: 

(A. 13) 1 - *(#043) = O(h) 

and 

(A- 14) *(£043) - tf(/*75o) = 0(h 2 ) 


Ro43 — 


a(uo — it 4) — /3hs- 


(A. 15) 
Similarly 
(A- 16) 

Therefore 

(A.17) 


= 1 


= 1 


( b - a)(u 3 — u 4 ) - (1 - (3)hs _ 

_ a(u 0 - u 4 ) + (6 - a)(u 3 - u 4 ) - 
(b - a)(u 3 - u 4 ) - (1 - /3)/is_ 

(au T + buy - s) — .bhbuyy — .5 Aau i:r — hbu xy — ,5/iSj 
(b — fl)uy — (1 — /?)s + 0(A) 


+ 0(h 2 ) 


= 1 - h- 


-j^bUyy .5 O/Uxx bu xy J, 


(6- a)u y - (1 - (3)s 


+ 0(A 2 ) 


^?750 = 1 — h- 


-.5bu yy + .5 au xx + .5 s x 
( b - a)u y - (1 - (3)s 

1 - #043 = 0(k) 


+ <D(h 2 ) 


and 

(A.i 8 ) 

Then because is Lipschitz continuous 

(A. 19) 1 -»(#043) = (1 -*(!)) + 0(A) 


R 750 ~ R 043 = h— jl — v + 0(A 2 ) = 0(h 2 ). 

(b — a)u y — (1 — (3)s v ' y ’ 
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and 

(A.20) V(R 75 o)-V(Ro43) = 0(h 2 ), 

which together with the fact that the 2D scheme is second order accurate proves the lemma. □ 


Appendix B. Fluctuation-splitting scheme for the inhomogeneous case 


Here we address the problem of maintaining linearity preserving property the NNL scheme on 
unstructured triangular meshes in the case of non-zero forcing term. Consider again the inhomo- 
geneous advection equation (4.1). The N scheme given by (4. 16), (4. 17) for the homogeneous case 
should be now modified. One target formula (Type I triangle, Fig. 7) now is 

(B.l) 5 3 ^ +1 = Szv% + At[ki{ui - u 3 ) + k 2 (u 2 - u 3 ) - a T ] 

and in two target case (Type II triangle, Fig. 7) 

= Sim" + - u 3 ) - 0a T ] 

= S 2 u% + At[k 2 (u 2 - u 3 ) - (1 - (3)<tt], 


(B.2) 


5j< +1 

S 2 u" +1 


where 

(B.3) 


0 = 


k i 

k\ + by 


and aj defined by (4.4). The following distribution formulae will give the inhomogeneous NNL 
scheme in the two target case 


(B.4) 


where 

(B.5) 


5iu" +1 = Siu* + At{[ k 1 (u i - u 3 ) - 0aj] 

+ ^(^132)[^2(^2 — W3) — (1 — P)&t]} 

5 2 u" +1 = S 2 «2 + At{[ k 2 {u 2 - u 3 ) - (1 - (3)ctt] 

~'HRl32)[k2(U2 ~ U 3 ) - (1 - 0)<?t]}, 


kijui ~ U 3 ) - 0OT 
k 2 (u 2 - u 3 ) - (1 - 0)<TT 


Identity 

(B.6) y(Ri32)[k 2 (u 2 - U 3 ) - (1 - 0)a T ] = ~ [Mm - « 3 ) - 0a T ] 

^132 

Positivity of this scheme can be demonstrated in the same way as for the homogeneous equation. 
The following result concerning the LP property of the constructed scheme can be formulated 


Theorem B.l. The NNL scheme (B.l), (B.2) is linearity preserving in the inhomogeneous case if 

(B.7) |*(i?)-l| <C|fl-l| 

for some positive constant C . 

Proof Define 

*«! = SiUj +1 - SjUj 

= At{[ki{Ui - U 3 ) - 0(Tt) + 'I'(#132)[M U 2 - « 3 ) - (1 - P)°t]} 

= P + V(R)Q 

(B.8) 

5u 2 = S 2 ti2 + 1 - S 2 U 2 

= At{[k 2 (u 2 - u 3 ) - (1 - 0)cr T \ - \P(f? 132 )[M “2 - u 3 ) - (1 - 0)a T ]} 

= Q-9(R)Q. 

The rest of the proof is the same as of Theorem 5.3. □ 
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To demonstrate the connection between the inhomogeneous NNL scheme and previously con- 
structed S’ scheme consider the grid illustrated in Fig.9. The fluctuations from the triangles A, B, C 
can be given by the following formulae 

4> A = { [1 - *(J2o 43)][(6 - o)(«3 - u 4 ) - (1 - P ) hs A ] 

+ a(no - U4 ) - @ hs ^ + ^(fio43)[(^> ~ a)(«3 - u 4 ) - (1 - 0 ) hs x ] } 

(B.9) ( t> B — — ajup - u ^) + (6 - a) (np - u 5 ) - /ls b ) 

4 > C = -|{ a ( u 7 - u 5 ) - 0 hs c + *(fl 75 o)[(& - a)(u 0 - « 5 ) - (1 - ( 3 ) hs c ] 

+ [1 - ^( j R 750 )][(6 - a)(rf 0 - m 5 ) - (1 - / 3 ) hs c ] } , 

where the underlined terms contribute to the residual at the node 0. Thus the update formula is 
“o +1 = u o + b( u 5 ~ u o) + a ( u 4 - u 5 ) 

-^\l/(f? 7 5o)[(&- a )( u 5 - u 0 ) + h ( 1 - /?)s c ] 

(/?043)[(h a )( u 4 - u 3 ) + h ( 1 - 0 ) s A ] 

(B.10) +^[s B + (l + 0)s A + (l-0)sc]}. 

Recalling the definition of s a ,sb,sc it is easy to see that (B.10) is identical to (A. 10). 
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